A Fast Algorithm for Generating Long Self-Affine Profiles 
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We introduce a fast algorithm for generating long self-affine profiles. The algorithm, which is 
based on the fast wavelet transform, is faster than the conventional Fourier filtering algorithm. In 
addition to increased performance for large systems, the algorithm, named the wavelet filtering 
algorithm, a priori gives rise to profiles for which the long-range correlation extends throughout the 
entire system independently of the length scale. 
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I. INTRODUCTION 



With the advent of the computer as a serious research 
tool, there has been a revolution in the quantitative de- 
scription of processes and structures that before were 
deemed too complex. Two of the key concepts used for 
this description are the fractal and its close relative, the 
self-affine structure [Q. In the early eighties, much ef- 
fort was spent identifying and describing various phys- 
ical systems having fractal or self-affine structure. As 
time went by, focus slowly shifted from pure description 
to asking why such structures would appear. This lead to 
the development of the science of complex growth phe- 
nomena. Now, many aspects of these are well under- 
stood. However, there are still hosts of interesting but 
unanswered questions lingering on — see e.g. Refs. (|-||] 
for recent reviews. More recently, focus has again be- 
gun to shift somewhat, and one sees work dealing with 
the physical consequences of the presence of fractal or 
self-affine structures. A concrete example of these three 
levels of development may be found in the study of frac- 
ture surfaces. In the early eighties, Mandelbrot et al. || 
characterized fracture surfaces as self-affine, in the early 
nineties attempts were made to understand why fracture 
surfaces are self-affine [||. Recently, phenomena such as 
two-phase flow in fracture joints have been studied 0. 

In order to study the physical consequences of the pres- 
ence of self-affine surfaces, algorithms generating these 
must be found. There are already several in existence, 
see e.g. Fedcr |IJ. However, subtle phenomena require 
the generation of huge surfaces. Two aspects of the al- 



gorithms then become important: (i) how self-affine are 
the surfaces that are generated, and (ii) how fast is the 
algorithm. 

The most popular algorithm used today is the Fourier 
filtering algorithm. This algorithm, which has a fast im- 
plementation thanks to the fast Fourier transform, con- 
sists of generating in the Fourier domain, uncorrelated 
Gaussian random numbers which are filtered by a decay- 
ing power-law filter of exponent — 2H — 1 where H is the 
Hurst exponent (to be defined in Section ||) . By taking 
advantage of the inverse fast Fourier transform, self-affine 
surfaces in real space, with the desirable correlations, are 
generated. 

It has previously been shown ||[)| that the Fourier fil- 
tering algorithm has the disadvantage that the self-affine 
correlations in the limit of large systems only exists over 
a fraction of the total system size. This is due to aliasing 
effects. For large enough systems this fraction might be 
well below 1%. One might overcome the above problem 
by e.g. temporarily generating a much larger surface than 
actually needed, and only use a small fraction of the total 
size. This is, however, not a very appealing approach, as 
the computer time and memory needed easily become too 
large. Another way of getting around this problem is due 
to Makse et al. [Q. Here the (Fourier-space) filter func- 
tion is modified by the introduction of a large momentum 
cut-off through the use of a modified Bessel function in 
the Fourier transform of the power spectrum. They show 
that this large momentum cut-off, while irrelevant for the 
large scale behavior in real space, is essential in order to 
suppress the aliasing effect and thereby obtaining sur- 
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faces with the desired scaling properties over the entire 
system size. 

In this paper we report on an alternative filtering al- 
gorithm based on wavelets which a priori, and without 
modifications, gives rise to self-affine correlations which 
extends (up to finite size effects) over the entire system 
independently of its size. This algorithm is also com- 
putationally cheaper than the traditional (or modified) 
Fourier filtering algorithm. 

This paper is organized as follows: In Section || we 
briefly review the defining properties of self-affine sur- 
faces. Here we also include some results which will prove 
useful to us later. Section III is devoted to the outline of 
the new algorithm. In Section [V we present numerical 



studies of the this algorithm. We conclude in Section 0. 



II. SELF-AFFINE SURFACES 

We limit our discussion in this paper to (1 + 1)- 
dimensional surfaces, which we will call profiles. A (sta- 
tistically) self-affine profile, h(x), is by definition a struc- 
ture which remains (statistically) invariant under the fol- 
lowing scaling relation 1 



x 
h 
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(2.1a) 
(2.1b) 



Here A is a real number and H, known as the roughness 
or Hurst exponent, characterizes this invariance. This 
exponent is usually in the range from zero to one. When 
H = 1/2, the profile is not correlated. An example of 
such a profile is the Brownian motion in one dimension. 
In this case, we interpret time as x and h(x) as the posi- 
tion on the Brownian particle at time x. When H > 1/2 
the profile is persistent, while when H < 1/2 it is anti 
persistent. 

We show in Fig. [I] an example of a self-affine profile 
generated by the algorithm to be presented in this pa- 
per. 

From the scaling relation (2.1), one can often relatively 
easily derive scaling relations for related quantities. In 
this paper we will later explicitly need the scaling rela- 
tion for the second order structure function 



S(Ax) = (\h(Ax + x) - h(x)\ 



(2.2) 



where (-) x represents the average over the position vari- 
able x, and the power spectrum P(q), defined as the 



Fourier transform of the height-height correlation func- 
tion. They scale as ]l]|| 



S(Ax) ~ {Ax) 2H , 



(2.3) 



and 



P(q) ~ q- 2H ~ l . (2.4) 
We will make use of these two scaling relations in Sec- 



tion 



IV 



III. THE WAVELET FILTERING ALGORITHM 

Recently, the wavelet transform jll]-[l3| has been used 
to analyze self-affine profiles 0,|l^l ■ The idea behind this 
analysis is as follows: We denote the wavelet transform 
of a function h{x) by W[/i](a, b), where a and b are the 
scaling and location variables respectively. They form 
together the wavelet domain. In Ref. [l^ the authors in- 
troduced what they called the average wavelet coefficient 
function, defined as IF[/i](a) = (|W[/i](a, b)\)b, where (■){, 
denotes the average over all the location parameters b 
corresponding to one and the same scale a. For a self- 
affine function, h{x), this quantity should scale as 



W[h](a) 



,H+l/2 



(3.1) 



In much the same way as the Fourier filtering algo- 
rithm is used for generating self-affine profiles via the 
fast Fourier transform, a wav elet based filtering tech- 
nique can be based on Eq. (3.1) in combination with the 
fast wavelet transform. The output of the fast wavelet 
transform [pj]-[T^] is a vector organized as a collection of 
various levels or hierarchies all of different lengths where 
each level, £, is associated with a corresponding scale 
ag. The two first components of this vector, also known 
as level £ = 0, are associated with the scaling function. 
All the other components are "true" wavelet coefficients, 
such that at level £, corresponding to scale at = 2~ £ , 
there are Nt — 2 l coefficients. These coefficients (using 
our convention) are arranged such that the coefficients of 
the highest level are found at the end of the vector, and 
the levels decrease monotonically towards the top of the 
vector, corresponding to the level I = 0. 

Hence, the wavelet based algorithm which we will be 
referring to as the wavelet filtering algorithm (WFA), 
consists of the following three steps: 



1 Strictly speaking, one should in order to fully define the 
self-affine structure, also introduce the topothesy which is de- 
fined as the length-scale I, over which the RMS-height, a, 
is a(l) = I. However, we will not explicitly need this latter 
quantity here, and will therefore simply neglect any further 
reference to it. 
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Generate in the wavelet-domain normalized un- 
correlated Gaussian numbers {iji}, with i = 
1,2, ... ,N where N is the number of discrete points 
Xi that together with hi = h(xi) constitute the self- 
affinc profile. N is assumed, due to the use of the 
fast wavelet transform, to be a power of 2. 



Filter these random numbers according to 



H+1/2 r\i 



(\v\h 



1,2, 



to obtain the wavelet coefficients {u>,}. Here am\ — 
2~ £ W represents the scale, at level £(i), where £(i) 
is defined as the level corresponding to the loca- 
tion index i of the vector Wi (or rji). Furthermore, 
(\v\)e(i) represents the average of the absolute value 
of those rji that together form level 

• Perform the inverse fast wavelet transform on {wi}, 
with the (compactly supported) wavelet of your 
choice, to obtain the (real-space) self-affine profile 
of predefined Hurst exponent H. 

With the wavelet filtering algorithm, good quality self- 
affine surfaces with predefined Hurst exponent can be 
generated. In Fig. [j] we show an example of a self-affine 
surface of Hurst exponent H = 0.6 and length N = 4096 
generated by the algorithm just outlined. It is worth not- 
ing that the above three steps can be modified in order to 
deal with surfaces in higher dimensions |l6| . In this case 
the speed of the surface generating algorithm becomes 
very important. 

One of the prominent features of the wavelet transform 
is that the basis functions, the wavelets, are localized in 
both space and frequency. This has as a consequence, 
among others, that there is no aliasing, or at least heavily 
suppressed as compared to the Fourier transform. This 
implies that the wavelet filtering algorithm should auto- 
matically result in surfaces which have the desired corre- 
lations over the entire length of the profile, and not just 
a small fraction of it. Hence, independent of the system 
size, the WFA is capable of generating profiles with well- 
defined long-range correlations. We will demonstrate the 
validity of this claim in Section [V. 



Before turning to the numerical studies of the WFA, 
we add some remarks regarding the computational effi- 
ciency of this algorithm. The most time-consuming part 
of the WFA is the inverse wavelet transform. To a good 
approximation, at least for larger system sizes, this time 
determines the overall computational time of the entire 
algorithm. The fast wavelet transform needs O(cN) op- 
erations, where c is a positive real number which value 
depends on the wavelet used (T^Jlj]]. Thus, the number 
of operations need for generating a surface by WFA is lin- 
ear in the number of points belonging to the profile. In 
comparison, the Fourier filtering algorithm, which speed 
is mainly controlled by the fast Fourier transform, needs 
0(iVlog 2 N) to generate a profile. For large system sizes 



the difference in execution time between WFA and the 
Fourier filtering algorithm becomes significant. 



IV. NUMERICAL RESULTS 

In order to test numerically the predictions of the pre- 
vious section, we have chosen to study the second order 
structure function, S(Ax), and the power spectrum P(q) 
of self-affine profiles generated with the wavelet filtering 
algorithm. The appropriate scali ng r elatio ns f or these 
two quantities are given by Eqs. (2.3) and (2.4). They 
will provide us with independent information enabling 
us to accurately quantify over which length scales the 
self-affine correlations exist. The numerical experiments, 
for which the results will be presented shortly, were per- 
formed as follows: We generated, by WFA, an ensemble 
of long self-affine profiles all with the same Hurst expo- 
nent H. For each profile the structure function and power 
spectrum were calculated, and these were averaged over 
the ensemble of profiles. 

In Fig. ^ we give the numerical results for the sec- 
ond order structure function obtained as described above. 
The predefined Hurst exponents, used by the surface gen- 
erator, were from bottom to top H = 0.8, 0.6, 0.4, and 
0.2 as indicated in the figure. The length of each pro- 
file was N = 2 25 = 33 554432. The number of profiles 
used in obtaining the averages was Nh = 50, and the 
wavelet used was of the Daubechies-type (-D12) [pTl|— p~3| . 
The dashed lines are regression fits to the numerical data. 
They corresponds (from bottom to top) to Hurst expo- 
nents of H = 0.80 ± 0.01, 0.60 ± 0.01, 0.41 ± 0.01, and 
0.22 ± 0.02, all consistent, within the errorbars, with the 
predefined exponents given above. One easily observes 
from Fig. || that the correlations extend over all scales 
except for the largest lags Ax. The reason that the last 
few large lags do not fit into this general picture is due 
to finite size effects. We have also undertaken the above 
analysis for different types of wavelets, taken from the 
Daubechies family, and for various system sizes, finding 
no results that are inconsistent with those presented in 

Fig. |. 

In order to make the comparison with the Fourier fil- 
tering algorithm even more apparent, we present in Fig. || 
the average power spectrum obtained using the same sur- 
faces as were used in Fig. |[ The correlations again span 
most scales. The dashed regression fits lead to the fol- 
lowing exponents (from bottom to top): H = 0.80±0.0f , 
0.61 ± 0.01, 0.41 ± 0.01 , and 0.20 ± 0.01, which again is 
in excellent agreement with the values of the Hurst ex- 
ponent used for the generation of the underlying profiles. 

Figs. U and || indicate that the self-affine correlations 
span all but the largest scales of the profiles. We stress 
that this is a generic property of the wavelet filtering algo- 
rithm, and no modification of the algorithm is needed in 
order to handle large system sizes in a satisfactory man- 
ner. This is a consequence of the celebrated property of 
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the wavelets being localized both in space and frequency. 

The calculations of this paper were performed on 
a SGI/Cray Origin 2000 supercomputer based on the 
R10000 chip from SGI. On this machine the average CPU 
time needed for generating a profile of the length used 
above (N = 2 25 ) was £\vfa = 45s and £ffm = 125s for 
the WFA and the traditional (or modified) Fourier fil- 
tering algorithm respectively. Hence the speedup gained 
by using the wavelet filtering algorithm over the Fourier 
filter algorithm is close to a factor 3. For system sizes 
N ~ 10 3 , we could not observe any significant different 
between the two algorithms. 



V. CONCLUSIONS 

To summarize, we have introduced a fast and simple al- 
gorithm for generating long (or short) self-affine profiles. 
This algorithm, named the wavelet filtering algorithm, 
is demonstrated to overcome the problem related to the 
aliasing effect which the traditional Fourier filtering algo- 
rithm is troubled with. Furthermore, the wavelet based 
filtering technique outperforms its Fourier-domain coun- 
terpart by large margins with respect to computational 
costs, at least for large system sizes. 
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FIG. 1. A self-affine profile, h(x), of length N = 4096 and Hurst exponent H = 0.6 generated by the wavelet filtering 
algorithm. The wavelet used in order to generate the profile was the D12 Daubechies wavelet. 

FIG. 2. The average second order structure function, S(Ax), obtained by averaging over N h = 50 samples (for given H) of 
the self-affine profile h(x) generated by WFA. All profiles were of length N — 2 25 = 33 5 5 4 432. The Hurst exponents used 
were (from bottom to top) H = 0.8, 0.6, 0.4, and 0.2, as indicated in the figure. The dashed lines are the best regression fits 
corresponding respectively to Hurst exponents (from bottom to top) H = 0.80 ± 0.01, 0.60 ± 0.01, 0.41 ± 0.01, and 0.22 ± 0.02. 



FIG. 3. The average power spectrum P(q) obtained by averaging over Nh = 50 samples (for given H) of the self-affine profile 
h(x) generated by WFA. All profiles were of length N = 2 25 = 33 5 5 4 4 32. The Hurst exponents used were (from bottom to top) 
H — 0.8, 0.6, 0.4, and 0.2, as indicated in the figure. The dashed lines are the best regression fits corresponding respectively 
to Hurst exponents (from bottom to top) H = 0.80 ± 0.01, 0.61 ± 0.01, 0.41 ± 0.01, and 0.20 ± 0.01. 
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